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Abstract 


The earth’s gravitational potential, as described by a spherical harmonic 
expansion to degree 180, has been compared to the potential implied by the 
topography and its isostatic compensation using five different hypothesis. 
Initially, series expressions for the Airy/Heiskanen topographic/isostatic model 
were developed to the third order in terms of (h/R) where h is equivalent 
rock topography and R is a mean earth radius. Using actual topographic 
developments for the earth, we found the second and third terms of the 
expansion contributed (on average) 30% and 3%, respectively, of the first term, 
of the expansion. With these new equations it is possible to compute depths 
(D) of compensation, by degree, using three different criteria: I) the power 

in the actual field at a specified degree should be the same as implied by the 
Airy/Heiskanen model; II) the topographic isostatic reduced potential should 
show minimum correlation with the earth’s topography; III) the norm of the 
residual potential should be a minimum. The results show that the average 
(over all degrees) depth implied by criterion I is 60 km while it is about 33 
km for criteria II and III with smaller compensation depths at the higher 
degrees. Another model examined was related to the Vening-Meinesz regional 
hypothesis implemented in the spectral domain. The fifth model tested took D 
to be a constant of 30 km at all degrees. We have compared these model 
fields with the actual field in terms of anomaly, geoid undulation and 
percentage differences, as well as with correlation coefficients and anomaly 
maps in the Caribbean/South America region. The differences between all 
models is small with the exception of the model defined by criterion I is used 
where larger differences are seen. For example, the average percentage 
differences between the OSU81 potential model and the five models outlined 
above, from degrees 15 through 180 is 87.2(1), 80.5(11), 79.8(111), 80.1(VM), 
80.5(D=30km). Finally oceanic and continental response functions are derived 
from the global data sets and comparisons made to locally determined values. 
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1. Introduction 


In recent years very efficient computer algorithms have been developed 
for spherical harmonic analysis up to high degrees and orders, (Colombo, 
1981). In addition, terrestrial gravity material has become much more 
complete, through the results of satellite altimetry, see e.g. (Rapp, 1981). 
These factB and the availability of a global set of l'xl* mean elevations make 
it very attractive to re-evaluate the older studies on global isostasy, carried 
out by Balmino et al. (1973) and others. A first step in this direction was an 
article by Rapp (1982). In Figure 1 of his paper a number of potential degree 
variance spectra are displayed. The spectrum of the "observed" gravity and 
the one derived from topographic heights and based on the Airy compensation 
model with a compensation depth of 50 km agree almost perfectly for the 
degrees 50 to 180. One may imply the following ideas from this figure: 

1. A compensation depth D=50 km seems to be much more realistic than the 
generally accepted one of D=30 km for the Airy model. 

2. Since the two spectra agree for n>50 but deviate for n<50 it may be 
feasible to apply an isostatic reduction to the observed gravity and, 
isolate those parts in the anomalous gravity field, that are not 
compensated isostatically, but are due to lateral inhomogenities in the 
earth's mantle. 

It is of general interest in geodesy, whether or not the potential of the 
isostatically compensated topography could be used as a device for the 
smoothing of the earth's gravity field, before it is used in geodetic 
approximation techniques, similar to the practice in local gravity field 
approximation (Forsberg and Tscherning, 1981). We realize that the full study 
of isostasy is a complicated issue. Isostatic behavior can vary from region to 
region and can be dependent on numerous factors that we do not consider. 
Nevertheless we feel that it is important to find out what we can learn from 
the models we formulate. 


2. Basic Equations for Airy Isostatic Model 


The gravitational potential at a point P outside the earth S is by Newton’s 
law of gravitation 


V(P) = G J d£ ( 


with G, the gravitational constant, p the density inside the earth, £ the 
distance between P and the infinitesimal volume element dS Q at Q. We insert 
for the inverse distance l/i PQ its spherical harmonic expansion: 

~T~ = IT 2 I?)" 9^77 ! Pnm (cosflp) P nm (cos0 Q ) [cosmAp cosmA Q 

* PQ r P n = o l "p J ^n+I m=0 

+ sinmAp sinmA Q ] (2) 


where P nm is the fully normalized associated Legendre function. 
In a shorter notation: 
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j“- = X fej" 2n+T Y nmo( (P) Y nroo< (Q) 

■* PQ r P n,m,o< lr V tn+l 


with 


Ynmot(P) = P(COS0 P ) 


cosmAp for o<=0 


( 3 ) 


sinmAp for <*=1 

If we further assume that the geocentric radius rp of P is greater than the 
radius of convergence of the series in (3), then (1) becomes after 
rearrangement 

V(P) = l rp("+0 Mfy J rg p(Q) Y nm0 <(Q) dS Q ) Y nfn0 <(P) 


= ^ (y" +1 (sdsii I (H" ' >(0) Y ""« (Q) ) Y -»« (p) (4 > 


with M, the mass of the earth and R its mean radius. From equation (4) 
follows the well-known expression 


GM 


y ( p ) = ^ I . 


&r 


' n moi 1 n md 


(p) 


with the dimensionless coefficients 

Cfin 


= c 


nm0( pR 3 (2n+l) 


hi ft)" p (Q) Y "««* (Q) d2 < 


(5) 


( 6 ) 


where M is replaced by 4/3npR 3 and p is the mean density of the earth. 

We consider now the model of local isostatic compensation by Airy . The 
coefficients of the potential generated by isostatically compensated topography 
Cjmd become, with equation (6): 

Cnrnc* = 2n+l) 4rr ~~ Y nmoi (Q) d<T q (7) 

The integration is performed over the unit sphere < r. The surface topography 
part is 

A T (Q) = [^j Pcr(Q) dr Q (8) 

with h, the surface elevation, and p cr the density of the topography in land 
areas. In ocean areaB p cr is replaced by p cr -p tu , where p m is the density of 
sea water. The compensation part in equation (7) is: 

A C (Q) = j R ° {%*) n+2 Ap( Q) dr q (9) 

J r=R-D-t J 

with D, the depth of compensation, t, the root (or anti-root) thickness, 
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Ap=p m -p cr > and p m the density of the mantle. If it is assumed that neither p crf 
or Pcr-Pu,, respectively, nor A p vary in radial direction, integration yields: 


*'«> * ^(0) & [m n+3 -l] 

(10) 

and 


AC (a , . M«, ^ [(¥)" +S - (S“)" + 1 

(11) 


In practical computations h(Q) is replaced by mean elevations usually l*xl*. 
For the geographical distribution of the densities one could in principle 
employ a classification into land, lake, sea, and ice areas. In practice, 
considering the limited quality of the data material, we can only distinguish 
between land and ocean areas with 


Pi = 


p cr = 2670 for hj>0 (land) 

Pcr-Pui = (2670-1030) for hj<0 (oceans) 


( 12 ) 


The derivation of the root thickness t shall be dealt with in the next chapter. 

With A T and A c evaluated e.g. in blocks of l*xl* and inserted into 
equation (7), the coefficients C£ m0 < of the isostatically compensated topography 
can be computed. Based on 5*x5* mean elevations such a computation was 
carried out by Lachapelle (1976) up to degree and order 36. 


3. Condition of Equilibrium 

The root thickness t(Q) is derived from the Airy model of local isostatic 
compensation. The equilibrium of mass condition between an element of mass 
at the surface and that at depth becomes for a spherical earth: 

fR+h fR-D 

p cr r 2 dr da - Apr 2 dr d<r 

J r=R J r=R-0-t 

or 


[(R+h) 3 - R 3 ] = ^ [(R-D) 3 - (R-D-t) 3 ] 

t 


We rewrite it as a third-order equation in 


R-D' 


fe)'-3fe)‘ + 3fe) --(£)' NT* 3(1)’ + 3(|)] =0 


or simply as: 

T 3 - 3T 2 + 3T - t,[H 3 + 3H 2 + 3H] = 0 
with T = jjip, H = jj, and tj = ^^(plp) • An iterative solution with the 
starting value T 0 = ijH yields up to the third order: 


(13) 


(14) 
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(15) 


T = To + AT = i,H[l + H(i,+1) - | H 2 (»j 3 -l)] 


or 


1 * % ( 515 )’ » I 1 + 5 ^ - 5 (D 2 *^- 1 )] (16) 

Since h/R < 2* 10“ 3 and »j 8 4.5 the error of a linear approximation remains 
below 1%. Hence we shall compute the root-thickness from 

t=£ ^(^)' h=4 - 46 ( ipy 2 ** <i7) 

with A P = Pn-Pcr = 3270 - 2670 £f. 

4. Series Expansion of Topographic and Isostatic Effects 

The problem with the direct application of formulas (10) and (11) for the 
computation of A T and A c , is, that they have to be computed anew at each 
block for each degree n. Hence, even with very efficient computer algorithms 
a high degree and order expansion, say n(max) = 100 or 200 becomes very 
tedious and time consuming. An alternative approach starts from a binomial 
expansion of the expressions in brackets of equations (10) and (11). Up to 
third order in h/R it is: 

*t(0) = „„«» jSj [l ♦ (n+3) wai „ In±32Inl21 (Ml) 1 + WWW 

(^r ♦ - l] 1 

* o.r(0) MO) [i + Sj 2 ^ + istmaill (M21)’] (is) 

A corresponding expansion of equation (11) yields 

*■<« ■ £ [(¥P - (¥)" +5 (i - <-) ^ 

+ (p» 3)(»+2) ||ffil)* . (m-3)(n+2)(n+l) + (()4) jj „ 

1 mo t«, m [i - ¥ 

The question is now at what term the series expansion has to be truncated 
without committing an unacceptably large error. 

We assume as an upper bound for h=10 km and take for R=6370 km, D=30 
km, and t * (p cr /Ap)h 8 4.45- h s 44.5 km. The resulting values of the first, 
second, and third order terms in the brackets of equations (18) and (19) are 
listed in Table 1 for a number of degrees n. Because of the chosen high 
value of h the numbers in Table 1 have the character of upper bounds. Still, 
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Table 1: Contribution of the second and third order terms in equations (18) 

and (19) relative to the first order term. 



A T Terms 

A c Terms 

n 

n+2 h 

(n+2) (n+1) h 2 

n+2 t 

(n+2) (n+1) t 2 

2 R 

6 R 2 

2 R-D 

6 (R-D) 2 

2 

0.003 

0.0 

0.014 

0.0 

10 

0.009 

0.0 

0.042 

0.001 

50 

0.041 

0.001 

0.186 

0.023 

100 

0.080 

0.004 

0.364 

0.088 

200 

0.158 

0.017 

0.721 

0.345 


in order to keep the truncation error below 10%, the second order terms have 
to be included for degrees greater than about 150 when computing A T , and for 
degrees greater than about 30, when computing A c . For A c one should even 
consider the inclusion of the third order for degrees higher than 150. The 
numbers in Table 1 are not very sensitive to changes in D. In conclusion for 
the high degree and order expansions, considered in this study the inclusion 
of at least the second order is required, in contrast to what is stated in 
(Dorman & Lewis, 1970, p. 3359). 

However, there is an additional aspect. We are actually interested in the 
combined effect of A T and A c , namely A T -A C , which enters into equation (7). 
If we insert for the root-thickness t from equation (17) 

t(8) - % (5§p h(Q) 

we find 

A'(Q) - A C (Q) = Pcr h {[l - Pjj£)"] + [i 

+ 1^!,.^ (¥r ]( !;n 

The alternating sign of the series expansion of 

(R-D1 n 

— g-J converges very slowly, it is true even 

for high degrees of n. The percentage values of the second, and third order 
terms in equation (20) relative to the first order are displayed in Figure 1. 
We see that throughout the entire spectrum from degree 2 to 200, the upper 
bound value (h=10 km) for the second order term either exceeds or remains 
little less than the first order term. The third order term reaches 25% of the 
first and second order for high degrees. If the upper bound value h=10 km 
is replaced by the approximate root mean square value of the topography of 


Pnr 

V 


(¥r] I 


( 20 ) 


A c causes the first order 
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about 2.5 km, the order of magnitude still remains between 18 and 46% of the 
first order, whereas the third order drops to below 1 % in this case. 

Now we insert equation (20) into equation (7), where ocean depths dj 
(which are negative in magnitude) are replaced by equivalent rock 
topography , (see ( 12 )): 


hj = P e r -£»u dj = 0.614 dj for hj < 0 ( 21 ) 

Then the (fully normalized) coefficients (up to and including 3rd order) are: 

= sir *T H 1 * Ml J ¥* (Q) d<r « + 

* ¥ [i ♦ X (¥)"”] «> d - 

♦ MM [ X - £ ( ¥ )-‘] A J him wa) *.} (22, 

With the surface spherical harmonic expansions 

h nmo( = 4 ^ J ^ Y nmc< (Q) d<r Q , 

<T 

h2„ m0( = ^ J ^r 1 Y nmc <(Q) d<r, and (23,a-c) 

<T 

h3 nmc( = £ J Y nmo( (Q) d<r , 


we obtain 


Cl 


nmo( 


- (¥ 11 » 


nmo< 


n+2 

2 


l 1 + ¥ (¥)”"1 


h 2 


nmc( 


+ ialiKa tn [i - pi (B^) ] h 3 „J . ( 24 ) 

This is the series expansion up to 3rd order for the computation of the 
potential coefficients of the isostatically reduced topography. In practice very 
often only the very convenient formulas of the first order expansion are used. 
They can be interpreted as a double layer expansion. But as we saw, this 
may result in considerable distortions in the computed coefficients. The need 
for a second order expansion was already emphasized in (Jung, 1950) with 
calculations carried out by Uotila (1965). In Arnold (1980) the second order 
expansion has been treated, too. The evaluation of the series expansion 
represented by (24) is much more efficient than the evaluation of the rigorous 
equations given by Lachapelle (1976). A comparison of coefficients computed 
from the rigorous application of ( 10 ) and ( 11 ) to the series results 
represented by (24) shows agreement on the order of 0.2xl0 -9 . 
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In numerical tests we noted that the power spectra of (h/R) 2 and (h/R) 3 
were approximately 10 -6 and 10 -13 , respectively, of the power of (h/R). 
However to see the whole picture we examined the relative contributions of the 
^nmai h2 nm0( , and h3 nm0( terms to the topographic-isostatic potential at various 
degrees. We found that the h2 nmc( contribution is on average 30% (a maximum 
of 110%) of the h nm0( contribution while the h3 nm0( terms contributes an average 
of 3% (maximum 11%) of the h nm0 ( contribution. These figures are consistent 
with the estimates made earlier. The significant role of the second terms is 
caused by the small value of the coefficient of h nin0 ( in (27). 

Another computation was the comparison of the potential implied by the 
three, two and one term topographic/isostatic models with D=30km. The 
anomaly, undulation, and percentage difference by degree were calculated. 
(See Section 8 for equations). The overal difference (degrees 2 to 180) 
between the three and one term models was 1.5m, 4.4mgal, and 35%. These 
differences were essentially the same (1.5m, 4.4mgal, 34%) for the two vs one 
term model. The magnitude of the undulation difference was greatest at the 
lower degrees while the anomaly and percentage differences were similar 
independent of degree. It is clear from these comparisons that a substantial 
error is caused by using only the first term in the model. For subsequent 
calculations, we will use the two term model. 


We finally note here that the series expression for the potential 
coefficients of the uncompensated topography can be written from (24) by 
letting D=R. We have: 

C nmo( = 2^1 ^ { h nmo( + ^ h2 nm0 ,+ I S+ ZhBtll h 3 nmo( } (25a) 

The coefficients of the isostatic compensation would be: 


qC _ O Pr.r 
nmo( 2n+l p 


f f H-D) "i-. 

p UrJ n ""“* 2 bp 

♦ is± ¥ a±u (¥)■"*“■ -} 


_ n+2 ^ | R-D j "~ 3 


h2 


nmao( 


(25b) 


The coefficients of the isostatically compensated topography would then be: 

nmol (25c) 


p* — p^ 

u nmo( = ' / nmo( ^ 


with equation (24) representing this difference. 


5. Preliminary Comparisons 

For the first series of numerical results we have computed the potential 
coefficient spectrum implied by the topography, and by the isostatically 
compensated topography. The spectrum of a set of coefficients A would be: 

*n (A) = E E C„L = E (C nm 2 + S n 2 ) (26) 

m <x m 

The elevations and depths for these computations have been taken from a 
l*xl* mean elevation file (TUG86.DTM) described by Siinkel (1986). This file is 
based on land elevations received from the Defense Mapping Agency Aerospace 
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Center in 1983 (with corrections for 12 values) and an averaging of a 5’x5’ 
data set (NGDC, 1983) over the oceans. In ocean areas the depths (d) were 
converted to rock equivalent topography using the nominal density of the 
crust and sea water using equation (21). No special consideration was given 
to ice areas because of lack of information. 

The harmonic coefficients of h/R and (h/R) 2 were computed to degree 180 
using program HARMIN (Colombo, 1981) using integrated associated Legendre 
functions. These coefficients are represented by equation (23, a,b,c). With 
these coefficients the potential coefficients of the uncompensated topography 
(Cj m0( ) and the isostatically compensated topography (C£ m0{ ) have been 
computed for the nominal compensation depth of D = 30 km. Plots of the 
spectra (from equation (26)) are shown in Figure 2. Also shown in this figure 
is the spectrum implied by the observed gravitational field of the earth as 
described by the OSU81 (Rapp, 1981) model which is complete to degree 180. 
From this figure we see the significant power in the uncompensated 
topography that is considerably reduced when isostatic compensation is taken 
into account. However the power of the D = 30 km case is less than that in 
the observed field which was noted by Rapp (1982). 

In Section 8 we will consider a number of comparisons between the 
observed potential field and model fields. At this point, however, we must 
discuss how the appropriate value of D is selected. 

6. Optimal Compensation Depths 

The isostatically reduced gravity potential dV is obtained by subtracting 
the potential due to the isostatically compensated topography, V 1 , from the 
"observed" gravity potential V: 

dV(P) = V(P) - Vi(P) (27) 

If the isostatic hypothesis of Airy fit reality perfectly, dV would be free of all 
influences of the earth’s crust and solely reflect the lateral density 

inhomogeneities of the earth's core and mantle. As mentioned in the 
introduction, from the degree variance spectra of V and of V 1 for the 

compensation depths D=30 km and 50 km in Rapp (1982), one might conclude 

that a compensation depth of 50 km is more realistic than one of 30 km. 

However, if one derives, from (27), the formula connecting the degree variance 
spectra of dV, V, and V 1 , we have: 

*S(V) = <r 2 (Vl) + 2cov n (V I , dV) + a 2 (dV) (28) 

One sees, that an agreement of <r 2 (V) and a^CV 1 ) does not necessarily imply 
that a„(dV) is small. One also has to take into account the covariance cov(V I , 
dV). 


We shall now discuss the computation of optimal compensation depths per 
degree , still within the Airy type model of compensation. Thus, the 
compensation depth D shall become a function of degree n. Three models for 
different optimality criteria shall be discussed. 
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Criterion 1: <rg(V) = gjlfV 1 ) 


As mentioned above this criterion has little physical relevance. Its application 
shall only result in an optimal agreement of the degree variances of the 
observed gravity field, <r„(V), and °f the topographic-isostatic model, ^(V 1 ). 
It does not necessarily produce a small residual field of dV. The criterion is 
included only for the purpose of showing that its application leads to 
compensation depths D n of about 50 km. To evaluate this definition we form 
the following function: 


f(D„) « <rg(V) - <r’(Vi(D n )) = 0 


with 

*S(V) = I l (C nm0 <)* 

m ot 

and with equation (24) 

«rS(Vl(D n )) = Z l (C^,*) 2 

m a 


(29) 

(30) 


= l l 


( JL Ecjl. 


m « '•2n+l p 


rU 1 - (¥“)>-« ♦ ^ 


] ^2 n m0 ( 


With A " = ShI 


+ iesptty _ £(K!) «]„,_}}■ 

= l l ( An[1 " Q " ]h " m0< + B "[ 1 + Q n" 3 ]h2 nm0( + 1' 

B n = A n 2±2, and Q n = BzLi 


(31) 


P 2 ’ “ R 

For the computation of the approximate compensation depths D„ the linear model 
for C* m0( is employed: 

^■n<no< ~ A n [l — Qn]h nfn0< (32) 

From 


f°(K) = ^(V) - Ap [1 - QSJ 2 I l h* m0( 

m & 

= <r*(V) - A’[l - QR]^(h) = 0 
which can be solved for D° in the following steps: 

! - Q n = /J EK 

" /A^(h) 


(33) 


m n - 




(h) 

‘ “I 1 - 

Based upon the approximate values D£, the actual depths can be obtained by 


(34) 
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Newton iteration. We write: 


D«+ l = D* - 


f(DJt) 

f'(D<) 


( 35 ) 


with the derivative f' of f from equations (30) and (31) and where the 
superscript i is an iteration index. 


f'(DA) 


3f I 
3D n Ida 


= 2 


1 1 

m & 


cA„,o<(dA){a„ | (^) n \ naM 


- B r 


£bjl 

Ap 


n-3 f R-DA 
R l R . 


n — 4 


and f' * 0. 



Using the OSU81 field and the rock equivalent topography expansions D n 
values have been computed. They are listed in Table 2, column 2 and 
displayed in Figure 3. The average D value is 64 km for degrees 2 to 180, 
and 58 km for degrees 30 to 180. 


Table 2. Depth (D) of Airy Compensation Based on Three Criteria 


Degree 

Criteria 

1 

2 

3 

Vening- 

Meinesz 

5 

144 km 

0 

0 

29 

10 

94 

65 

64 

29 

15 

54 

40 

40 

29 

20 

40 

27 

27 

30 

25 

40 

31 

31 

30 

30 

34 

22 

23 

30 

50 

52 

30 

31 

31 

70 

49 

29 

30 

31 

90 

60 

30 

30 

32 

110 

61 

35 

37 

32 

130 

79 

38 

40 

33 

150 

56 

29 

32 

34 

170 

68 

25 

29 

34 

175 

61 

31 

33 

35 

180 

56 

26 

29 

35 


Criterion 2: lcov„(V T ,dV) I = min 

A very valid requirement is the choice of the depths of compensation such 
that the topographic-isostatic reduced potential dV shows minimum correlation 
with the earth's topography. The correlation can be expressed as: 

cov n (V T , dV) = I I Crj mc( dC nm0 ( (37) 

m at 
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DEGREE 

Figure 3. Compensation Depths for the Three Criteria Considered as a Function 
of Spherical Harmonic Degree. 


or using equations (25c): 

cov n (V T , dV) = l l Cj m0( (C 

nmo( ^nmof ^nmof) 

m oi 

= cov n (V, VT) - <r„(V T ) + cov n (V T , V c ) (38) 

In linear approximation equation (37) becomes with the definitions in (31): 

cov n (VT, dV(Dp) ) = A n cov n (V, h) - A*<r*(h) + A’Q"<rji(h) (39) 

If by definition D„ = {D„l 0 < DJ < R} then A^Q n <r^(h) is a monotonously 
decreasing function with 


ASQ-^(h) = { AS S s<h > f" jjs : l 


R 

Now three cases can be distinguished in order to attain the minimum of 
cov n (VT, dV(D®)): 

if A n cov n (V, h) - A„<r„ (h) >0 => D° = R 
if A n cov n (V, h) - Ap<r„ (h) < A*<r*(h) => D° = 0 
otherwise, lcov n (V T , dV(D„)l = min implies: 
cov n (V T , dV(D°)) = 0 


(40) 
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Then, with equation (39): 

Q „ _ Ag<rg(h)-A n cov n (V, h) _ cov n (V, h) 
y A^(h) A n <r„(h) 

or 

ds = «[t - (i - x:»?h ) h) )' / 1 (4i) 

For the actual solution, again the three cases D n = R, D n = 0, and 0 < D n < R 
can be distinguished. In the latter case, i.e. 0 < D n < R, the compensation 

depths are derived from a Newton iteration, as before, with f(D n ) = 
cov n (V T , dV) according to equation (38). 

The resulting compensation depths are listed in Table 2, column 3, and 
displayed graphically in Figure 3. The average D value is 32 km for degrees 
2 to 180, and 31 km for degrees 30 to 180. 

Criterion 3: «dV» 2 - min 


Whereas criterion 2 seems to be more relevant for geophysicists, the next 
criterion (3) is especially interesting for geodetic application. In geodesy 
topographic-isostatic reduction is mainly applied to produce a smooth residual 
field which allows simple and accurate interpolation. This is achieved with 
criterion 3. It is: 


B dVH 2 = X {l X [C nm0( - CiU + Cfj m o((Dn)p} 

n d J 

= X (x X [C-Ua + (CT m0( ) 2 + (Cnmo((D n ) ) 2 - 2 C nmo{ CT m0< 

n v m o( 

+ 2^nmc^nmo<(I)n ) “ 2C J m£ *C fj mo < ( D n ) ] j . 

HdVH 2 = min, is obtained from a solution of 

l l ai j“ + °nmo< ^nmof J U 

or X X dC nmo( = 0 

m d au r\ 

3C c 

In linear approximation, — with equation (25b), can be written: 

®u_ 



A n § Q n_1 h nm0 ( 


(42) 


(43) 


- - 2 Qn~ l Q T 

pj ^ mo( 

Consequently equation (43) becomes: 


- | Q n-1 cov n ( V T , dV) = 0 
Equation (44) is satisfied either by 


(44) 
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I: Q"-i = (*p)"“ = o 

which is true for D® = R, or by 
II: cov n (V T , dV) = 0 

which is identical to the linear case of criterion 2, treated above. Solution I 
leads to C{; m0 < = 0, and results in a pure Bouguer reduction, compare equation 
(25). In contrast to solution II, solution I represents a maximum of equation 
(42) and is therefore of no interest. Hence, in linear approximation, criteria 2 
and 3 agree, or in other words the minimum norm of dV and minimum 
correlation of V T and dV result in one and the same compensation depths D„. 

In second and higher order approximation, iidV» 2 = min. differs from 
lcov(V T , dV)l = min. The solution is again obtained by Newton iteration with 
equation (42) in the following form: 

f(D - ) = 5 5 dC -”“ 

= • « *" 8 8n + 4 ^ 2nmo ). 

{c nm o( - A n (l-Q")h nmt * - B„[l + Q n— 3 J h2 nm0 ( + — J 

f(D„) = Q n-1 {- | A n [cov n (V, h) - A n (l-Q")<r 2 (h) - B n 

(! + Q n-3) CO v n (h, h 2 )] + 2^ ^ BnQ -3 

[cov„(V, h 2 ) - A n (l-Q n )cov n (h, h 2 ) ~ B n (l + Q"” 3 ) <r 2 (h 2 )) 

(45) 

from which f'(D n ) = -rr— follows by straightforward differentiation. 

dU n 

The resulting compensation depths are given in Table 2, column 4 and 
displayed in Figure 3. They differ from those obtained from criterion 2 
typically by 1 to 2 km with the greatest differences at the higher degrees. 
The average D value is 33 km for both degrees 2 to 180 and 30 to 180. 

At this point we have three different procedures to compute the depth of 
compensation by degree. Given the depths, which will depend to some extent 
on a potential and topographic model of the earth, the topographic/isostatic 
potential can be computed using equation (24) with D now a function of n. 
Numerical comparisons between these and other models will be discussed in 
Section 8. 

7. Optimal Vening Meinesz Isostatic Models 

Half a century ago Vening Meinesz argued that the simple isostatic concept 
of Airy/Heiskanen, which is governed by a strictly local compensation, is not 
very realistic because, due to its elasticity, the mantle is able to support a 
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certain amount of local topographic load without local yielding, whereas in the 
Airy/Heiskanen model free mobility between vertical mass columns is 
presupposed - a highly unlikely assumption (Vening Meinesz, 1939). 
Therefore, he proposed a regional compensation model which is now denoted 
the Vening Meinesz model (Heiskanen and Moritz, 1967). 

Regional mass compensation can be described by a smoothing operator, 
applied to the non-smoothed isostatic masses implied by the Airy/Heiskanen 
model. This idea has been pursued by Sunkel (1985, 1986). In the following 
we shall describe how to best estimate both the depth of compensation and the 
parameter (s) of the smoothing operator. 

In linear approximation the harmonic coefficients of the topographic/ 
isostatic potential for the Airy/Heiskanen model are given by 

C " m = 2n+l £ p C t 1 ~ (~) ] h " m <* ( 46 ) 

(cf. eq. (24)), where the term with "1" is the contribution of the topography, 
and the term with ((R-D)/R) n is the contribution of its isostatic compensation. 

If the isostatic masses are subject to a smoothing procedure, implied by 
a homogeneous and isotropic smoothing operator B, 

B(P,Q) = Z (2n+l) /?„ P n (cos^p Q ) (47a) 

n=o 

with eigenvalues 0 n , 

= 2n P B(t) P n (t)dt (47b) 

J -i 

(P n ... Legendre polynomial; t := coat//] if/ ... spherical distance), we then obtain, 
by observing the convolution theorem, the harmonic coefficients of the 
topographic/isostatic potential in linear approximation for a Vening Meinesz 
model, 

C " m = 2n+l P ~p t 1 _ (“r”) ^"3 hnmc( * ( 48 ) 


(Note that the topographic part remains unchanged; only the isostatic part 
changes because of smoothing!) 

Various models for the set of eigenvalues {/9 n } could be envisioned, each 
particular model controlling the features of the corresponding smoothing 
operator 0. We list here a very few: 

a) /?„ = 1 V n 

In this particular case there is no smoothing involved. Therefore, this 
operator represents the standard Airy/Heiskanen model with the 
compensation depth D as the only parameter. 

b) /Jo = 1 J n = 0 V n > 0 

In this case the smoothing operator degenerates into the global average 
operator with equal weight; the smoothed compensation masses form a 
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homogeneous layer which is equivalent to a point mass at the origin. 
Therefore, this operator represents the other extreme case of the 
Airy/Heiskanen model where the "level" of compensation is at the origin, 
D = R. 

c) = 0 V n 

This set of eigenvalues represents the annihilation operator which 
annihilates all isostatic masses. Therefore, only topographic masses are 
left (unchanged) in this particular case. Consequently, the corresponding 
harmonic coefficients are those of the topography only, which represent 
the topographic potential. 

d) With 



a smoothing operator is implied which is identical to the isostatic model 
used in section 6: a model where the compensation depth formally depends 
on the degree n, rather than being constant as in the general case of eq. 
(48). We think it is important to stress and therefore, we repeat it: the 
multi— compensation depth model discussed in section 6 represents, in 
linear approximation, also some type of smoothing operator which is 

applied to the compensation masses. 

It should be obvious that there exists an infinite number of other models 

for the eigenvalues {/?„} of the smoothing operator B, e.g. e -b n or the like; 
(here b is the only model parameter of the smoothing operator). Generally, 
the smoothing operator of eq. (47a) depends on a certain number p of model 
parameters; together with the compensation depth parameter D we have 

therefore to estimate p+1 parameters of the isostatic model. 

The three optimization criteria for the estimation of the model parameters, 
which have been discussed in section 6, can be applied here as well. In 

principle the estimation procedures agree with those of section 6 with two 
small differences: 

a) In the case of the Vening Meinesz model each harmonic C£ m0 < depends on 

all p+1 model parameters, while in section 6 a harmonic of degree n 

depended exclusively on the degree-specific compensation depth D n : 

b) through the smoothing procedure, the operator B enters into higher 

degree terms like h2 nmc( , h3 nm0< , ... in a rather complicated manner, 

requiring for its computation a (quickly converging) iteration procedure. 

Because of property (a) the estimation of the degree-specific compensation 
depths D n according to criteria I - III of section 6 is possible on a degree by 
degree basis. Here, in the general smoothing model approach, this is no 
longer possible (exception: according to model (d)). Instead of the 

solution of N max (= highest degree of model) optimization problems with one 
parameter D n , we have to solve now one optimization problem with p+1 
parameters. 
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Because of property (b) an iteration procedure is involved: eq. (48) 

presupposes a linear relation between the harmonic coefficients of the 
topographic-isostatic potential and the topographic height. But according to 
eq. (25b) this relation is non-linear because of the higher degree terms h2 nm0( , 
h3 nm0( , ... (Note that in the Vening Meinesz model h nm0( is replaced by /9 n h nm0( ; 
therefore h J nm0{ is the harmonic coefficient of the smoothed topography B*h, 
raised to the power J). Therefore, the following iteration process offers itself: 

1 = 0 . 
get a priori estimates for the isostatic model parameters Dl 0 ' and 

b(°) (here b is the vector comprising all p parameters of the chosen 
smoothing operator); harmonic analysis of the topography and its low 
powers, yielding h nm0 <, h2 nm0( , ...; 
i := i + 1 

smoothing of compensation masses, implied by B( ,_1 )*h using DV'“ l ) 
and b(’” 1 ), yielding a smoothed topography h(0; 

harmonic analysis of low powers of h( 1 ) yielding a first order term 
Oji') of eq. (48) which will be considered as a function of D and b, 
and higher order (correction) terms 02^*), ••• which will be 

considered as constants in the subsequent step; 

optimal estimation solution (least-squares adjustment or least-squares 
collocation with parameters) for D and b with Taylor point 
b0“ l ) and linearization restricted a to the^ linear term C^ 1 ) of step 2; 
stop if IDt 1 ) - d( ,-1 )i < eq and lb(0 - b( ,-1 )l < E b , else go to (1). 

The result of this iteration process yields both best estimates for the 
isostatic model parameters (compensation depth D and the parameters of the 
smoothing operator B), and the set of harmonic coefficients of the topograhic - 
isostatic potential of the implied Vening Meinesz model. 

This procedure has been applied to the determination of a two parameter 

model with D = constant and = e -1 * 2 " 2 . The result is D=29 km and b=0. 00223. 

The values of D n found by setting 0 n = e -62 " 2 to ( (R-D n )/(R-D) ) n are shown in 
Table 2. The values are similar to the other values above degree 20. 


step 0: 


step 1: 
step 2: 


step 3: 


8. Observed and Model Comparisons 

The two main results obtained so far are: 1. In isostatic computations 
linear approximation may result in significant distortions. Second and higher 
order approximations can be implemented very easily. 2. In addition to the 
classical Airy/Heiskanen local compensation model with constant compensation 
depth D, three models with degree dependent compensation depths D n and a 
Vening-Meinesz regional compensation model have been derived. 


We now turn to a comparison of the different models. In Figure 4 we 
show the spectra of the OSU81 (Rapp, 1981) field, the model with degree 
dependent D values based on criterion 3 and the Vening-Meinesz model 

with smoothing coefficients /9 n = e -t,2n2 (D = 29 km, b = 0.00223). The D=30 km 
model, if plotted, would be almost the same as the Vening-Meinesz model. All 
models, except the one based on criterion 1, have less power than the observed 
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Figure 4. Log of the Potential Degree Variances of the OSU81 Potential 
Coefficient Solution, the Potential Implied by the 
Topographic/Isostatic Model Using Criterion 3, and the Potential 
Implied by the Vening-Meinesz Model, as a Function of Spherical 
Harmonic Degree. 


potential at the higher degrees. 

Next a set of quantities is introduced that can be used to compare two 
sets of potential coefficients. We define the difference between the 
coefficients C nmoi of the observed gravity potential and the isostatic 
coefficients C£ m0 ( (model) as derived from the mean elevations by applying one 
of the five models as 

ACnmoi ~ ^nmd — C nmc< (model). (49) 


The first quantity of interest is the root mean square undulation differences 
between degrees n x and n 2 : 

<5N = [r* 2 2 I l AC* m0( ]^ (50) 

L n=ni m « ■» 

The next is the root mean square anomaly difference between degrees n x and 
n 2 : 

6g = [y 2 I* (n-1) 2 l l AC 2 **]* (51) 

L n-n x m & J 

We next have the percentage difference by degree: 


P„ 


l AC 2 *,* 

*JS(V) J 


% 

100 


The average percentage difference between degree n x and n 2 is: 
1 


P = 


n 2 

l p. 


n 2 -ni+l n =n j n 


(52) 


(53) 


The ffn(V) are the dimensionless signal variances of the observed gravity 

potential as derived from the OSU81 set of spherical harmonic coefficients. As 
can be seen from the above formulas at each degree the relative behavior 
among the AC nm0 < derived from the five models remains the same in eqs. (50), 
(51), and (52). Only when the root mean square (r.m.s.) quantities in a range 

from degree n, to n 2 are computed this relative behavior could change. This 

is due to the multiplication by (n-1) 2 in (51) or the division by <r„(V) in (52), 
which changes the relative weight of certain degrees or degree ranges, 

respectively. 

The root mean square undulation and anomaly differences are chosen, 
because they provide a good insight into the expe.cted size of the residual 
undulation and anomaly field. The percentage difference is a valuable measure 
if the differences are small as compared to ffp(V). Its disadvantage is that it 
is unbounded, i.e. the percentage can reach infinity. 


A parameter containing the same information as the percentage difference 
is the smoothing coefficient s„. It has been used by Tscherning (1985) in 
comparing different geopotential models to the potential implied by the 
topography. It is defined as 
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s 


n 


Z Z 

= m U 

ff n(V) 


(54) 


Another measure is the correlation (coherence) coefficient by degree 
between C nmot and Cj mc< (model) 


£ E ^nmQ( ^nmo( (model ) 
<r„(V)<r n (V I (model)) 


(55) 


Finally we have the average correlation coefficient between degree n, and n 2 


P = 


Z 2 P„ 


n 2 -n x +l n=ni 


(56) 


The degree correlation coefficient (or coherence coefficient) is a quantity that 
can be used to judge the agreement between any two quantities represented 
in a spherical harmonic expansion. Our particular case can be explained by 
taking the linear approximations to C^ m0 < from equation (24): 

= di *T I 1 - (¥)1 h "«« 

so that 

- 2 ^*? ( l -'(¥!"] “" <h) 

Hence, in linear approximation the correlation coefficient per degree becomes, 
with (54): 

Z Z C nm0( Cj mC ( (model) Z C nmo(^nmc( 

a - (^7') 

Pn <r n (V)<r n (VI (model)) <r n (V)<r n (h) 

The latter is the correlation coefficient between the observed gravity potential 
and the observed equivalent rock topography. We conclude that in linear 
approximation, the />„ are independent of the choice of the isostatic model. 
Consequently the differences in the correlation coefficients of the five models 
must be small and represent the influence of the second and higher order 
terms. In addition, since the second order term of C£ m0 < (cf. (24)) is positive 
by definition, one can show that the correlation coefficients of the linear 
models must be higher than those for the corresponding quantities in higher 
order approximations. On the other hand, since we know that the models in 
linear approximation represent rather unrealistic double layer models, we 
conclude that higher correlation coefficients do not necessarily imply a more 
realistic model. 


Table 3 shows the anomaly, undulation and percentage differences, and the 
correlation coefficients for the various topographical isostatic models. The 
simple Airy model is represented through a D = 30 km depth of compensation. 
Other depths (28 and 32 km) were used but no significant changes noted. In 
addition, three degree ranges (2 to 180, 15 to 180, 30 to 180) have been used 
for the comparisons. All models yield very similar comparisons except for the 
case when D is computed by Criterion 1. This points out the clear problem of 
seeking D values that yield the same power as in the observed field. 
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Table 3. Comparisons Between 0SU81 Model and Various Topographic/Isostatic 
Models for Three Different Degree Ranges. 


Solution 

Anomaly Diff 
(mgals) 

Undulation Diff 
(meters) 

Percentage 

Difference 

Correlation 

Coefficient 

2 

180 

15 

180 

30 

180 

2 

180 

15 

180 

30 

180 

2 

180 

15 

180 

30 

180 

2 

180 

15 

180 

30 

180 

D=30 

D=Crtl 

D=Crt2 

D=Crt3 

Vening- 

Meinesz 

20.64 

24.11 

20.16 

20.15 

20.58 

16.27 

17.64 

16.27 

16.12 

16.20 

15.43 

16.74 

15.42 

15.30 

15.36 

30.57 

45.35 

29.59 

29.59 

30.57 

2.30 
2.48 

2.31 
2.26 

2.29 

1.46 
1.58 

1.47 
1.44 

1.46 

81.48 
89.19 

80.48 
80.43 

81.12 

80.51 

87.18 

80.54 

79.75 

80.12 

80.41 

87.18 

80.38 

79.75 

80.00 

.575 

.593 

.584 

.587 

.581 

.599 

.616 

.603 

.603 

.604 

.599 

.617 

.602 

.603 

.605 


Computations were also made with the TUG87 (Wiser, 1987) elevation model for 
selected cases. We found that this model yields better agreement with the 
0SU81 field than the TUG86 model. For example, using TUG87 the average 

percentage difference for D=30 km, degree 15 to 180, is 79.47% compared to 

80.51% with TUG86. The correlation coefficient, in the same case is 0.614 as 
compared to 0.599. All computations were not carried out with TUG87 as it 
became available after most of the computations were completed for this paper. 

Table 4 contains information on the smoothness coefficient (s n ) for selected 
degrees and for the same three degree ranges used previously. The smaller 
the s n value, the better the agreement between the observed field and the 
models. All isostatic models, except that based on criterion 1, yield basically 
the same smoothing coefficients with those from criteria 2 and 3 being slightly 
better, i.e. smaller, for the 2 to 180 degree range. It lies in the very 
definition of criterion 3 - minimum norm of dV = V - V 1 - that the smallest 

values of AC nm< * for each degree and therefore also of fiN, <5g, P n , and s n must 

be obtained with this criterion. In practice the differences between criteria 2 
and 3 and the Vening-Meinesz model are very small. Obviously the smallest 
covariance values cov n (V T , dV) must be obtained for criterion 2. We note that 
the use of the TUG87 model would have yielded slightly smaller P n values. 
For example, for degree range 30 to 180, the s value is 0.626 for TUG87 vs 
0.647 for TUG86. 

Correlation coefficient information by degree range has been shown in 
Table 3. The correlation coefficients by degree show no significant variation 
between the models as expected. At degree 12 the correlation is about 0.04 
which is unusual because most correlations in this degree area are on the 
order of 0.6. This small correlation at degree 12 also holds for other potential 
coefficient models. 

Anomaly degree variances are a measure of the spectrum in the anomaly 
domain. We write: 
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Table 4. Relative Mean Square Coefficient Differences (s) Between 0SU81 and 
Various Topographic/Isostatic Models 



Fixed D (km) 

Variable D 

Vening 

Meinesz 

t 

30 

1 

2 

3 

2 

1.05 

2.86 

.98 

.98 

1.06 

5 

1.18 

2.78 

.97 

.97 

1.18 

10 

.68 

.63 

.53 

.53 

.68 

12 

1.18 

1.89 

1.05 

1.05 

1.16 

15 

.55 

.58 

.50 

.50 

.55 

20 

.76 

.86 

.75 

.75 

.75 

25 

.46 

.51 

.46 

.46 

.46 

30 

.73 

.79 

.67 

.67 

.72 

50 

.64 

.78 

.64 

.64 

.64 

70 

.59 

.73 

.59 

.59 

.59 

90 

.70 

.90 

.70 

.70 

.70 

110 

.60 

.67 

.58 

.58 

.58 

130 

.65 

.76 

.63 

.63 

.64 

150 

.68 

.79 

.68 

.68 

.67 

170 

.79 

.95 

.79 

.79 

.78 

180 

.67 

.79 

.67 

.66 

.66 

2-»180 

.664 

.796 

.648 

.647 

.658 

15->180 

.648 

.760 

.649 

.636 

.642 

30-»180 

.647 

.760 

.646 

.636 

.640 


c n = r 2 (n-l) 2 l I C„ m0 ( (58) 

ffl of 

The c n values computed from the models used in this paper are interpreted to 
refer to a sphere corresponding to the mean earth radius of 6371 km due to 
the spherical approximations of the model. Values of c n are given in Table 5. 
The model with the most power is the variable D case with criterion 1. The 
other models have power by degree, and cumulatively, that is significantly 
less than what is in the observed field as represented by the 0SU81 field. 

The coefficients have next been used to create gravity anomaly maps in 
the Carribean region and in northern South America. This area was selected 
due to the complex nature of the gravity field due to the presence of both a 
trench area (the Puerto Rican Trench) and a substantial mountain range (the 
Andes Mountains). From the point of view of an area where isostatic 
compensation would be most apparent, this area may not be ideal. These maps 
have been created using the potential coefficients from degree 30 to 180 of 
the 0SU81 field (Figure 5) and the Vening-Meinesz regional compensation field 
(Figure 6). The anomaly plots for the D=30 km case and the depths from 
criterion 1 case are very similar. 

The OSU81 field reflects the known anomaly field structure in this region. 
There is clearly a correlation of the significant anomaly structure between the 
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0SU81 field and the topographic/isostatic model. However the range of the 
anomalies is much smaller in the model case. 

Table 5. Anomaly Degree Variances (mgal 2 ) 




Fixed D 

Variable D 

Vening- 

Meinesz 

t 

0SU81 

30 km 

1 

2 

— 

3 


15 

3.0 

0.9 

3.0 

1.1 

1.6 

0.8 

30 

2.5 

2.0 

2.4 

1.1 

1.1 

1.9 

50 

3.6 

1.4 

3.6 

1.4 

1.5 

1.4 

70 

2.7 

1.1 

2.6 

1.1 

1.1 

1.2 

90 

2.6 

0.9 

2.6 

0.8 

0.9 

0.9 

100 

2.6 

0.8 

2.6 

0.9 

0.9 

0.8 

120 

2.1 

0.8 

2.6 

1.1 

1.3 

0.9 

140 

2.2 

0.8 

2.2 

1.0 

1.0 

0.8 

160 

1.8 

0.8 

1.8 

0.8 

0.9 

0.9 

180 

1.6 

0.8 

1.5 

0.7 

0.7 

0.9 

2->180 

586 

176 

582 

201 

210 

181 

15->180 

412 

169 

408 

173 

182 

174 

30-»180 

372 

149 

368 

155 

164 

155 


Quantities that are primarily used in the geophysical literature are the 
continental and oceanic response functions. The response function (as 
discussed here) is a measure of the transfer of information between 
topography and gravity field quantities , such as Bouguer or free air gravity 
anomalies. Discussions of response functions in flat earth approximation 
applying two dimensional Fourier transforms can be found, for example in 
(Dorman & Lewis, 1970; McKenzie & Bowin, 1976, or McNutt, 1979). A general 
review is given by Boekelo (1985). A discussion of the response function 
technique for a spherical earth can be found in (Dorman & Lewis, ibid), 
(Boekelo, ibid), (Cazenave et al., 1986) and other sources. 

The response function technique assumes a linear relationship between 
gravity anomalies and equivalent rock topography. Let us define the 
coefficients of the free air gravity anomaly field as 

£>nmo< - 7(n — l)^nmo( (59) 

and those of the Bouguer gravity anomaly field (compare (25a)) as 
£nmo( = y(n— 1) {^n»o( — ^nmo(} 

= r(n-l) {c nm0( - 2^1 ^ h nm0( ) (60) 

Then the assumption is that there exists an unknown linear transfer function 
Z n (g f ,h) such that 
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LATITUDE 





LATITUDE 



Figure 6. Free-Air Gravity Anomalies Implied by the Vening-Meinesz Model 
from Degree 30 to Degree 180 in Northern South America and the 
Caribbean Sea; Contour Interval = 20 mgal. 
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6nmc< - Z n (g f »h) h nmc( + eJjMji 


(61) 


and a transfer function Q„ (g b ,h) such that 


£nmo( ~ Qn(^ b > ^)^nmo( + e nmoi 


(62) 


In both cases the e£ m0( and e b m0( are assumed to be small and (almost) 
randomly distributed. Then the coefficients Z n and Q n can be determined in 
the least-squares sense from 


Z n (g f , h) 


and 


Z Z Snmof^n md 


I I hnmofknmd 

m 6t 


Q n (g b » h) 


Z Z S b mo(^nmo( 

- m u 

Z Z hnmdhnmd 

m d 


(63) 


(64) 


In function Z n is usually referred to as oceanic or free-air response, whereas 
Q n is denoted isostatic or continental or Bouguer response . 


and 


The theoretical values of Z n and Q n for our five models become 

z » * s?i *f I 1 - (¥)"'■] 

0" = r(n-l) ^ ^ (Sz2)Y 


where 


(}„ = 1 for the Airy/Heiskanen model (D = const) 
fR-D 1 n 

f) n = I— ■ - j r a | for the models based on criteria 1, 2, and 3 and 


(65) 

( 66 ) 


2 2 

P n - e -b n for the Vening Meinesz model. 

Figure 7 shows the oceanic (free-air) response functions Z n of our five 
models and as directly derived from the data. It can be seen that they all 
agree very well, except the one based on criterion 1. As to be expected, the 
free air anomaly field has little correlation with the topography at low 
frequencies whereas for high frequencies the correlation increases. It is 
remarkable that the response as derived from the data stays significantly 
below that determined from the models for wavelengths smaller than 300 km. 


In Figure 8 the analogous curves of the continental (Bouguer) response 
functions are given. Again the curve based on criterion 1 departs from all 
others. For the long wavelengths down to 2000 km all other curves show 
approximately a Q n value of -0.10 mgal/m, somewhat lower than the well known 
Bouguer gradient of -0.1119 mgal/m. For small wavelength the correlation of 
the Bouguer anomalies with topography decreases, as to be expected. 
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Figure 7. The Oceanic (Free-Air) Response Function, Z„ (in 0.1 mgal/m) By 
Degree, for Various Topographic/Isostatic Models, and as Derived 
From the 0SU81 Model. 
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Figure 8. The Continental (Bouguer) Response Function Q n (in 0.1 mgal/m), by 
Degree, for Various Topographic/Isostatic Models and as Derived 
From the OSU81 Model, with Sample Values From Regional Studies. 
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Also shown on Figure 8 are sample values of continental response of 
Eastern U.S., Western U.S., and Australia, as derived from regional studies. 
The values are taken from (McNutt, 1978 and 1980). The values are scattered 
around the curves of the global models and data. This would indicate that 

our data can be seen as a global average. Only for wavelengths smaller than 
about 400 km are the sample values throughout smaller, which could indicate 
that especially for these wavelengths there exists a significant difference in 
response between continental and oceanic lithosphere. It should be pointed 
out again, that a general drawback of the response technique is, that it 
assumes a linear relationship and consequently implies a simple, double layer 
approximation. 

9. Conclusions 


This paper has considered improved procedures to compute the 
topographic/isostatic effects due to Airy compensation taking into account 
three terms in a series expansion. Normally only the first term is used. 
Tests showed that the second term in an expansion contributes an average of 
30% of the first term while the third term contributes an average of 3%. 
Overall differences between one and three component models of ± 4.4 mgal and 
35% indicate the importance of the higher order terms. The potential field 
implied by this theory has been compared with an observed field defined by a 
set of potential coefficients to degree 180. Comparisons were also made with a 
Vening-Meinesz type model derived by Sunkel. Comparisons were carried out 
for fixed depths of compensation and for depths computed as a function of 
spherical harmonic degree using three different criteria. Two of these criteria 
(II and III) led to nearly the same depths of compensation while criterion I 
led to depths of compensation larger than normally accepted for the Airy 
hypothesis. The results of the comparison of the theoretical models with the 
observed field were given in Tables 3, 4, 5. No substantial differences occur 
between the models except for the case of depth Criterion I which implies 
substantially poorer agreement with the observed field. 

The rather high compensation depths D n , for the very low degrees, that 
were obtained, when applying criteria 1, 2, and 3 (see Table 2) could be an 
indication that the compensation of the large scale topographic features cannot 
be explained by some simple Airy or Vening-Meinesz type of model. The great 
depths could be a hint that the long wavelengths compensation is related to 
the convection mechanism in the upper mantle. 


The Vening-Meinesz model applied by Sunkel results in the response 
function (cf. (54) and (65)): 


q - y ( n -i) — — £j tr fO] % 

l) 2n+l Up I R J P " 


(67) 


with p n = e“ b n . Turcotte & Schubert (1982) discuss the flexure of the 
lithosphere under periodic loading, derived from the equation of equilibrium of 
a deformed plate. The result - when transformed from a plane to a spherical 
approximation - is a smoothing factor 
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where D' is the flexural rigidity (see also e.g. (Banks et al., 1977)). Hence 
this smoothing model behaves basically as l/(l+kn 4 ) with k a constant. It can 

easily be shown that a smoothing model with e -b4n4 and an appropriate choice 
of b, can be made to agree almost perfectly with the above expression. This 
allows us to interpret the constant b of the Vening-Meinesz regional 
compensation model in terms of the flexural rigidity D' of the the deformed 
plate model in (Turcotte & Schubert, ibid). 

Four of the five models discussed here show a very similar overall 
behavior, despite the different assumptions on which they are based. This 
leads to a twofold conclusion. First, a rather large variation in the chosen 
isostatic models has little effect on the resulting isostatic gravity potential. 
Second, the other way round, despite today’s availability of rather good 
observed gravity and topographic elevations, the data are not selective 
enough to identify one particular isostatic model. This also implies that the 
uncertainty of parameters derived from the isostatic models, when dealing with 
global models. It means also that one of the two primary objectives, stated in 
the introduction can hardly be met, namely to isolate and display those parts 
in the anomalous gravity field, that are not compensated isostatically, but are 
due to lateral inhomogeneities in the earth’s mantle. This basically reconfirms 
Dahlen (1982; p. 3947) who states "... that this can never be done 
unequivocally even if the topography d and density anomalies p i are known 
exactly, since it requires that the concept of local isostasy be defined 
extremely precisely, more precisely, in fact, then the concept itself probably 
warrants." 


The remaining gravity potential, after subtracting the isostatic part from 
the observed one, must to some extent be attributed to the fact, that the 
approach chosen here does not easily allow the application of different 
isostatic models for different tectonic zones (e.g. continents and ocean areas) 
and that no corrections were made for geoid-age, depth-age effects or 
sediment loading as has been done e.g. in (Cazenave et al., 1986) or (McNutt 
and Shure, 1986). 

Since the isostatic behavior of the earth is dependent on a number of 
factors, and considering that such behavior varies substantially from area to 
area, global models cannot be expected to reflect the full picture. This paper, 
however, has examined what can be learned from global modeling techniques, 
and does not imply that the mechanisms employed are the only ones involved 
in the complex physical process of isostasy. 

Our results and conclusions depend to a lesser extent on the accuracy of 
the topographic model and the observed potential field model. The 
topographic model used is a clear improvement over previous models. However 
we still do not have ice thickness used in the computations. The observed 
potential field could be updated to more current models (e.g. Rapp and Cruz, 
1986, Wenzel, 1985) but no substantial changes are expected. 
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